COMMENT ⊗   VALID 00005 PAGES
C REC  PAGE   DESCRIPTION
C00001 00001
C00002 00002	_What to do Until the Doctor Comes_
C00004 00003	Numerical Precision
C00007 00004		I	3 digit SUM	8 digit SUM	Exact SUM
C00012 00005	Numerical Precision Exercise
C00014 ENDMK
C⊗;
_What to do Until the Doctor Comes_

We are talking about a doctor of philosophy specializing in numerical analysis.
In everyday life, everyone should know a little about medicine - not so that
they can prescribe antibiotics for themselves, but so that they can recognize
the symptoms that need to be checked out by a professional.

Just as our bodies are inhabited by micro-organisms that normally are harmless,
but occasionally are lethal, our programs are full of micro-errors arising from
the limited precision of computer arithmetic.  Occasionally these small errors
during a computation lead to a catastrophically large error in the result.  This
is the professional domain of the numerical analyst.  Like the proctologist's,
his job is not to everyone's liking.

We are not going to try to make junior numerical analysts out of you, but we
are going to show you what problems to look for; a competent programmer must
be aware of these dangers.  According to Prof. Wm. Kahan, several major
passenger aircraft are plagued with aerodynamic problems arising from faulty
numerical analysis in the design process.

++++++++
Write More.
Numerical Precision

Computers normally work with numbers represented by a limited number of digits.
In the early days of electronic computing, many of them worked with decimal
numbers; typically, integers were ten digit numbers, and real numbers were the
eight leftmost nonzero digits, together with two digits of coded information
about decimal point location.  Now most computers use binary (base - 2) 
arithmetic, but the consequences of using limited precision are qualitatively
much the same.  Let's look at a simple example, a program that forms the first N
multiples of 1/3, as it would run on a decimal computer.

	SUM:=0;
	FOR I:=1 TO N DO
		SUM:=SUM+1/3

If the computer only kept the three most significant digits of every real number
(unrealistic, but it keeps the arithmetic simple), the calculated value of
1/3 would be 0.333, and the variable SUM would become 0.333, 0.666, and 0.999.
At the next addition, the true sum 1.332 would be truncated to 1.33, and so
long as SUM was under 10, it would increase at each step by 0.33.  When SUM
reaches 10, it starts increasing by 0.3, and when it reaches 100, it doesn't
change at all.

The table below sums up the history of the computation using such three digit
real numbers, using eight digit real numbers, and using the exact numbers.
	I	3 digit SUM	8 digit SUM	Exact SUM

	1	   0.333	0.33333333	  0.333...
	3	   0.999	0.99999999	  1.0
	4	   1.33		1.33333333	  1.333...
       30	   9.91		9.99999991	 10.0
       31	  10.2	       10.333332  	 10.333...
      300	  90.9	       99.999909        100.0
      301	  91.2	      100.33323	        100.333...
      330	  99.9	      109.99980		110.333...
      331	 100.	      110.33313		110.333...
		(no further change)
     3000        "	      999.99090	       1000.0
     3001	 "	     1000.3242	       1000.333...
    30002	 "	     9999.7575	      10000.666...
    30003	 "	    10000.090	      10001.0
   300273	 "	   100000.00	     100091.0
  3027545	 "	   999999.76	    1009181.666...
  3027546  	 "        1000000.0	    1009182.0
 33027546	 100	 10000000.	   11009182.0
			 (no further change)

As you can see, the three digit sum, when it reaches 100, is wrong by about
10%, even though the error committed at each step is only about 0.1%; after
that the sum becomes stuck, and its value is useless. Similar phenomena happen
eventually even with eight digit computation.  After about 30000 steps, the value
of the sum is only correct to four places.  After some 33 million steps, only
the first digit of the sum is correct, and after that the value sticks.

While 33 million steps is a vast computation for people, for present computers
it can be over in a minute or less.  It is important, then, to watch for situations
where the limited, though high, precision of the computer gives a very imprecise
result.

Usually, such imprecise results happen in one of two ways: accumulation of small
errors, or effects of cancellation.  In cumulative computations like the one we
have seen, a rough rule of thumb is that the cumulative error is the precision
of the computer (typically 10↑{-8}) times the number of steps in the computation;
if a computation accumulates a sum of more than 10000 addends, one should be
wary of trusting its accuracy.

The other form of dangerous computation entails subtracting two nearly equal
numbers to get a much smaller number.  The result of such a step may be very
imprecise, and may contaminate the rest of the computation.

********************
(Need example here of a cancellation error)

To estimate the effect of errors in a computation, study the algorithm as
though it were being carried out by malicious employees who try to make the
results inaccurate without being noticed.  In order not to be noticed, at each
step the result of that step can differ from a correctly calculated value by,
say, one part in 10↑{8}.  Each result is in effect multiplied by some number
1+E,|E|<10↑{-8}, with a different  E  at each step.

If I evaluate  A+(B*C)-D,  say, what the malcontents give me is 
(B*C)(1+E↓1)instead of B*C; (A+(B*C)(1+E↓1))(1+E↓2) instead of A+B*C; and
((A+(B*C)(1+E↓1))(1+E↓2)-D(1+E↓3).  

Expanding this and ignoring products of E's, we get
		
	(A+B*C-D)+E↓1(B*C)+E↓2(A+B*C)+E↓3(A+B*C-D).

If B*C or A+B*C is very much larger than A+B*C-D, as can happen with
cancellation, the error may be a large fraction of the result even though the
E's are small.
Numerical Precision Exercise

The evaluation of

1-1/2+1/3-1/4+....+1/9999-1/10000

was carried out by these four program fragments, with the results shown

beside them.  The correct result is

S:=0.0;
FOR I:=1 TO 5000 DO
	S:=(S+1/(2*I-1))-1/(2*I)

S:=0.0;
FOR I:= 5000 DOWNTO 1 TO
	S:=(S+1/(2*I-1))-1/(2*I)

SPLUS :=0.0;SMINUS:=0.0;
FOR I:=1 TO 5000 DO
	BEGIN
	SPLUS:=SPLUS + 1/(2*I-1);
	SMINUS:=SMINUS + 1/2*I)
	END;
S:=SPLUS-SMINUS

SPLUS:=0.0; SMINUS:=0.0;
FOR I:=5000 DOWNTO 1 DO
	BEGIN
	SPLUS:=SPLUS + 1/(2*I-1);
	SMINUS:=SMINUS + 1/(2*I)
	END;
S:=SPLUS-SMINUS

The calculation was done in decimal arithmetic, keeping the six most significant
digits of each number (that is, numbers were truncated, not rounded, to six
digits).

The errors differ enormously.  Explain why the errors are the approximate size
they are.  Notice that the first program is the obvious one to write, and
has much poorer precision than the others.